#######################################################################################################
## Replication Script                                                                                ##    
## State Capacity, Refugee Reception, and Attitudes towards Refugees                                 ## 
## Study 2: Expanding Reception Capacity Improves Attitudes towards Ukrainian Refugees               ##  
#######################################################################################################

#######################################################################################################
#### packages ####
library(rio)
library(ggplot2)
library(gridExtra)
library(ggpubr)
library(patchwork)
library(scales)
library(car)
library(dplyr)
library(tidyr)
library(udpipe)
library(stm)
library(Rtsne)
library(rsvd)
library(geometry)
library(modelsummary)
library(marginaleffects)
library(sandwich)
library(clubSandwich)
library(lmtest)


set.seed(0987)

#######################################################################################################

#######################################################################################################
#### sessionInfo() ####

# R version 4.4.0 (2024-04-24)
# Platform: aarch64-apple-darwin20
# Running under: macOS Sonoma 14.5

# Matrix products: default
# BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib  
# LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

#locale:
#  [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

# time zone: America/Chicago
# tzcode source: internal

# attached base packages:
#  [1] stats     graphics  grDevices utils     datasets  methods   base   

# other attached packages:
# [1] lmtest_0.9-40          zoo_1.8-14             clubSandwich_0.5.11    sandwich_3.1-1         marginaleffects_0.21.0
# [6] modelsummary_2.1.1     geometry_0.5.2         rsvd_1.0.5             Rtsne_0.17             stm_1.3.7             
# [11] udpipe_0.8.11          tidyr_1.3.1            dplyr_1.1.4            car_3.1-2              carData_3.0-5         
# [16] scales_1.4.0           patchwork_1.2.0        ggpubr_0.6.0           gridExtra_2.3          ggplot2_3.5.2         
# [21] rio_1.2.0           

# loaded via a namespace (and not attached):
# [1] tidyselect_1.2.1    farver_2.1.2        R.utils_2.12.3      fastmap_1.2.0       bayestestR_0.15.2   digest_0.6.36      
# [7] estimability_1.5.1  lifecycle_1.0.4     NLP_0.3-2           magrittr_2.0.3      compiler_4.4.0      rlang_1.1.6        
# [13] tools_4.4.0         data.table_1.16.0   knitr_1.48          ggsignif_0.6.4      labeling_0.4.3      stopwords_2.3  
# [19] xml2_1.3.6          RColorBrewer_1.1-3  abind_1.4-5         tinytable_0.3.0     withr_3.0.2         purrr_1.0.2      
# [25] R.oo_1.26.0         grid_4.4.0          datawizard_1.0.0    fansi_1.0.6         xtable_1.8-4        tm_0.7-16      
# [31] future_1.33.2       globals_0.16.3      emmeans_1.10.3      insight_1.0.2       cli_3.6.5           mvtnorm_1.2-5  
# [37] generics_0.1.3      rstudioapi_0.16.0   future.apply_1.11.2 performance_0.12.2  parameters_0.24.1   magic_1.6-1    
# [43] stringr_1.5.1       parallel_4.4.0      effectsize_0.8.9    matrixStats_1.3.0   vctrs_0.6.5         Matrix_1.7-0      # [49] slam_0.1-55         rstatix_0.7.2       listenv_0.9.1       glue_1.8.0          parallelly_1.37.1   codetools_0.2-20 
# [55] cowplot_1.1.3       stringi_1.8.4       gtable_0.3.6        tables_0.9.28       tibble_3.3.0        pillar_1.11.0    
# [61] htmltools_0.5.8.1   R6_2.6.1            lattice_0.22-6      R.methodsS3_1.8.2   backports_1.5.0     SnowballC_0.7.1  
# [67] broom_1.0.6         Rcpp_1.0.13         coda_0.19-4.1       checkmate_2.3.1     xfun_0.46           pkgconfig_2.0.3 
#######################################################################################################

#######################################################################################################
#### Figure 4 Data Import ####

gt_data <- rio::import("Hlatky_UKR_Ref_Study2_GoogleTrend_Data.csv")

gt_data$Day <- as.Date(gt_data$Day, format = "%m/%d/%y")


specific_date <- as.Date("2022-05-03")



gt_plot <- ggplot(gt_data, aes(x = Day, y = Value, linetype = Term)) + geom_line() + scale_x_date(date_labels = "%d %b", date_breaks ="2 days") + theme_minimal() + geom_vline(xintercept = as.numeric(specific_date), linetype = "dashed") + xlab("Date") + ylab("Standardized Search Frequency") + theme(legend.position="bottom") + theme(legend.title=element_blank()) + ylim(50,100)

gt_plot <- gt_plot + guides(linetype = guide_legend(nrow = 2))

gt_plot

# ggsave("Fig4_GT.jpeg", gt_plot, width = 6.75, height = 5, dpi=600)

#######################################################################################################

#######################################################################################################
#### STM Data Cleaning and Preprocessing ####

text_data <- rio::import("Hlatky_UKR_Ref_Study2_Text_Data.csv")

ud_model <- udpipe_load_model("slovak-snk-ud-2.5-191206.udpipe")

annotations <- udpipe_annotate(ud_model, x = text_data$text)

anno_df <- as.data.frame(annotations)

anno_df$lemma <- ifelse(
  grepl("^[A-Z0-9\\-]+$", anno_df$token),
  anno_df$token,
  anno_df$lemma
)

content_words <- subset(anno_df, upos %in% c("NOUN", "ADJ", "VERB","NUM","ADV","PROPN"))

cleaned_texts <- content_words %>%
  group_by(doc_id) %>%
  summarise(text = paste(lemma, collapse = " ")) %>%
  ungroup()

sk_stopwords <- stopwords::stopwords(language = "sk", source = "stopwords-iso")

processed <- stm::textProcessor(documents = cleaned_texts$text,
                                metadata = cleaned_texts,
                                language = "",  
                                customstopwords = sk_stopwords)

out <- prepDocuments(processed$documents,
                     processed$vocab,
                     processed$meta,
                     lower.thresh = 2)

#######################################################################################################

#######################################################################################################
#### Study 2: Figure A6 STM Diagnostic Plot ####

storage <- searchK(documents = out$documents,
                   vocab = out$vocab,
                   data = out$meta,
                   K = 2:8,
                   seed = 1234,
                   init.type = "Spectral")
graphics.off()

# pdf("FigA6_STM_Diagnostic.pdf", width = 7, height = 8)


plot(storage)

#dev.off() 

#######################################################################################################

#######################################################################################################
#### Study 2: Table A6 STM Topic Top Words with English Translations ####

stm_model <- stm(documents = out$documents,
                 vocab = out$vocab,
                 data = out$meta,
                 K = 4, 
                 init.type = "Spectral",
                 seed = 1234)

semcoh <- semanticCoherence(stm_model, out$documents)
exclus <- exclusivity(stm_model)

topic_prevalence <- colMeans(stm_model$theta)
names(topic_prevalence) <- paste("Topic", 1:length(topic_prevalence))
print(topic_prevalence)

summary(stm_model)
#######################################################################################################

#######################################################################################################
#### Study 2: Table A7: Study 2 Descriptive Statistics ####

data <- rio::import("Hlatky_UKR_Ref_Study2_Survey_Data.csv")

data$datum <- as.Date(data$datum, format = "%m/%d/%y")

# recode so higher values greater problem perception, greater acceptance, more trust 
data$q16_H_problem_Ukrainian_refugees <- car::recode(data$q16_H_problem_Ukrainian_refugees, "1=4; 2=3; 3=2; 4=1")
data$q10_P_accept_Serb_Romanians <- car::recode(data$q10_P_accept_Serb_Romanians, "1=4; 2=3; 3=2; 4=1")
data$q4_R_trust_standard_media <- car::recode(data$q4_R_trust_standard_media, "1=4; 2=3; 3=2; 4=1")

data$treatment_full <- as.factor(ifelse(data$datum<"2022-05-03",0,1))
table(data$treatment_full)

datasummary(q16_H_problem_Ukrainian_refugees + q4_R_trust_standard_media + q10_B_who_power_does_not_matter + q10_S_want_anti_corruption + q10_P_accept_Serb_Romanians + q17_A_beneficial_sending_troops_nato + q18_D_feel_informed_EU + D2pov_age + r11_income +  y1_left_right+ y2_lib_con + c3_religious_attendance + D5_size_settlement + as.factor(treatment_full) + as.factor(D1_gender) + as.factor(D3_education) + as.factor(D6_region)  ~ Mean + SD + Min + Max + N, data=data, dinm=FALSE, output = "latex")

#######################################################################################################

#######################################################################################################
#### Study 2: Table 3: Pre- and Post-Treatment Balance ####

means <- data %>% group_by(treatment_full) %>% dplyr::summarize(
  D1_gender = mean(D1_gender),
  D2pov_age = mean(D2pov_age),
  D3_education = mean(D3_education),
  c3_religious_attendance = mean(c3_religious_attendance),
  r11_income = mean(r11_income, na.rm=TRUE),
  y1_left_right = mean(y1_left_right),
  y2_lib_con = mean(y2_lib_con),
  D5_size_settlement = mean(D5_size_settlement),
  D6_region= mean(D6_region))

pre <- data[data$treatment_full==0, ]
post <- data[data$treatment_full==1, ]

t.test(pre$D2pov_age, post$D2pov_age)
t.test(pre$D1_gender, post$D1_gender)
t.test(pre$D3_education, post$D3_education)
t.test(pre$r11_income, post$r11_income)
t.test(pre$y1_left_right, post$y1_left_right)
t.test(pre$y2_lib_con, post$y2_lib_con)
t.test(pre$c3_religious_attendance, post$c3_religious_attendance)
t.test(pre$D6_region, post$D6_region)
t.test(pre$D5_size_settlement, post$D5_size_settlement)


#######################################################################################################
#### Study 2: Figure 5: Pre- and Post-Treatment Avg. Values of UKR. Refugee Problem Perception ####

aggregated_data <- data %>% group_by(datum) %>% dplyr::summarize(avg = mean(q16_H_problem_Ukrainian_refugees), na.rm=TRUE)

aggregated_data$treated <- ifelse(aggregated_data$datum<"2022-05-03", 0,1)

specific_date <- as.Date("2022-05-03")

timeseries_plot <- ggplot(aggregated_data, aes(x = datum , y = avg, fill=as.factor(treated))) + geom_point(shape = 21, size = 3, alpha=20) + scale_x_date(date_labels = "%d %b", date_breaks ="2 days") + geom_vline(xintercept = as.numeric(specific_date), linetype = "dashed")+theme_minimal() + scale_fill_manual(labels=c("Untreated", "Treated"), values = c("white", "lightgray")) +  theme(legend.position="bottom") + theme(legend.title=element_blank()) + xlab("Date") + ylab("UKR. Refugees are a Problem for SK")

timeseries_plot

# ggsave("Fig5_TimeSeries_UKR.jpeg", width = 6.5, height = 5, dpi=600)

#######################################################################################################

#######################################################################################################
#### Study 2: Figure 6 Effect of EU Announcement on UKR. Refugee Problem Perception and Placebo Checks ####

outcome_model_bivariate <- lm(q16_H_problem_Ukrainian_refugees ~ as.factor(treatment_full), data=data)
summary(outcome_model_bivariate)

outcome_model_control_imbalanced <- lm(q16_H_problem_Ukrainian_refugees ~ as.factor(treatment_full) + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con , data=data)
summary(outcome_model_control_imbalanced)

outcome_model_control_all <- lm(q16_H_problem_Ukrainian_refugees ~ as.factor(treatment_full) + D2pov_age + as.factor(D1_gender) + as.factor(D3_education) + r11_income + y1_left_right + y2_lib_con + c3_religious_attendance + as.factor(D6_region) + D5_size_settlement, data=data)
summary(outcome_model_control_all)

data$treatment_full <- as.factor(ifelse(data$datum >= "2022-05-01" & data$datum < "2022-05-03", 0,
                                        ifelse(data$datum >= "2022-05-03" & data$datum < "2022-05-05", 1, NA)))

range2_model <- lm(q16_H_problem_Ukrainian_refugees ~ as.factor(treatment_full) + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=data)
summary(range2_model)

data$treatment_full <- as.factor(ifelse(data$datum >= "2022-04-30" & data$datum < "2022-05-03", 0,
                                        ifelse(data$datum >= "2022-05-03" & data$datum < "2022-05-06", 1, NA)))

range3_model <- lm(q16_H_problem_Ukrainian_refugees ~ as.factor(treatment_full) + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=data)
summary(range3_model)

data$treatment_full <- as.factor(ifelse(data$datum >= "2022-04-29" & data$datum < "2022-05-03", 0,
                                        ifelse(data$datum >= "2022-05-03" & data$datum < "2022-05-07", 1, NA)))

range4_model <- lm(q16_H_problem_Ukrainian_refugees ~ as.factor(treatment_full) + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=data)
summary(range4_model)

alt_data <- data[data$datum!="2022-05-03", ]

alt_data$treatment_full <- as.factor(ifelse(alt_data$datum<"2022-05-04",0,1))


drop_may3_model <- lm(q16_H_problem_Ukrainian_refugees ~ as.factor(treatment_full) + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=alt_data)
summary(drop_may3_model)

alt_data_2 <- data %>% filter(!datum %in% c("2022-05-02", "2022-05-03", "2022-05-04"))

alt_data_2$treatment_full <- as.factor(ifelse(alt_data_2$datum<"2022-05-05",0,1))


drop_may2_4_model <- lm(q16_H_problem_Ukrainian_refugees ~ as.factor(treatment_full) + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=alt_data_2)

summary(drop_may2_4_model)

alt_data_3 <- data %>% filter(!datum %in% c("2022-05-01", "2022-05-03", "2022-05-05"))

alt_data_3$treatment_full <- as.factor(ifelse(alt_data_3$datum<"2022-05-06",0,1))


drop_may1_5_model <- lm(q16_H_problem_Ukrainian_refugees ~ as.factor(treatment_full) + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=alt_data_3)
summary(drop_may1_5_model)

vcov_list_itt_models <- list(
  vcovHC(outcome_model_bivariate, type = "HC2"),
  vcovHC(outcome_model_control_imbalanced, type = "HC2"),
  vcovHC(outcome_model_control_all, type = "HC2"),
  vcovHC(range2_model, type = "HC2"),
  vcovHC(range3_model, type = "HC2"),
  vcovHC(range4_model, type = "HC2"), 
  vcovHC(drop_may3_model, type = "HC2"), 
  vcovHC(drop_may2_4_model, type = "HC2"), 
  vcovHC(drop_may1_5_model, type = "HC2")
)

itt_models_labelled <- list("Drop May 1-5 Obs."=drop_may1_5_model,"Drop May 2-4 Obs."=drop_may2_4_model, "Drop May 3 Obs."=drop_may3_model, "4-Day Bandwidth"=range4_model, "3-Day Bandwidth"=range3_model,"2-Day Bandwidth"=range2_model, "All Controls"=outcome_model_control_all, "Imbalanced Controls"=outcome_model_control_imbalanced, "Bivariate"=outcome_model_bivariate)

model_plot <- modelplot(itt_models_labelled, size=.1,  vcov=vcov_list_itt_models,  coef_omit = c(-2), conf_level=.95, color="black", facet = TRUE, coef_rename = "", ) + geom_vline(xintercept = 0, linetype = 'dashed') + theme_minimal() + xlim(-0.3,0.3) + xlab("Treatment Effect") +  theme(axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 8), axis.title.x = element_text(size = 8.5)) + labs(title="ITT Estimates + Robustness") + theme(plot.title = element_text(size = 7, hjust = 0.5)) + theme(plot.title = element_text(size = 10))

data$treatment_full <- as.factor(ifelse(data$datum<"2022-05-03",0,1))

placebo_outcome_model_1 <- lm(q10_P_accept_Serb_Romanians ~ as.factor(treatment_full) + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=data)
summary(placebo_outcome_model_1)

placebo_outcome_model_2 <- lm(q10_S_want_anti_corruption~ as.factor(treatment_full) + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=data)
summary(placebo_outcome_model_2)

placebo_outcome_model_3 <- lm(q17_A_beneficial_sending_troops_nato ~ as.factor(treatment_full) + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=data)
summary(placebo_outcome_model_3)

placebo_date_data <- data[data$datum < as.Date("2022-05-03"), ]

placebo_date_data$treatment_full <- ifelse(placebo_date_data$datum<"2022-04-28",0,1)

placebo_date_model_1 <- lm(q16_H_problem_Ukrainian_refugees ~  as.factor(treatment_full) + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=placebo_date_data)
summary(placebo_date_model_1)

placebo_date_data$treatment_full <- ifelse(placebo_date_data$datum<"2022-04-29",0,1)

placebo_date_model_2 <- lm(q16_H_problem_Ukrainian_refugees ~  as.factor(treatment_full) + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=placebo_date_data)
summary(placebo_date_model_2)

placebo_date_data$treatment_full <- ifelse(placebo_date_data$datum<"2022-04-30",0,1)

placebo_date_model_3 <- lm(q16_H_problem_Ukrainian_refugees ~  as.factor(treatment_full) + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=placebo_date_data)
summary(placebo_date_model_3)

placebo_date_data$treatment_full <- ifelse(placebo_date_data$datum<"2022-05-01",0,1)

placebo_date_model_4 <- lm(q16_H_problem_Ukrainian_refugees ~  as.factor(treatment_full) + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=placebo_date_data)
summary(placebo_date_model_4)

placebo_date_data$treatment_full <- ifelse(placebo_date_data$datum<"2022-05-02",0,1)

placebo_date_model_5 <- lm(q16_H_problem_Ukrainian_refugees ~  as.factor(treatment_full) + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=placebo_date_data)
summary(placebo_date_model_5)

vcov_list_placebo_models <- list(
  vcovHC(placebo_outcome_model_1, type = "HC2"),
  vcovHC(placebo_outcome_model_2, type = "HC2"),
  vcovHC(placebo_outcome_model_3, type = "HC2"),
  vcovHC(placebo_date_model_1, type = "HC2"),
  vcovHC(placebo_date_model_2, type = "HC2"), 
  vcovHC(placebo_date_model_3, type = "HC2"), 
  vcovHC(placebo_date_model_4, type = "HC2"), 
  vcovHC(placebo_date_model_5, type = "HC2")
)

placebo_models_labelled <- list("DV: Serbs/Romanians"=placebo_outcome_model_1, "DV: Corruption Eval."=placebo_outcome_model_2, "DV: NATO Involvement"=placebo_outcome_model_3, "04/28 Placebo" =placebo_date_model_1, "04/29 Placebo"=placebo_date_model_2, "04/30 Placebo" =placebo_date_model_3,"05/01 Placebo" =placebo_date_model_4, "05/02 Placebo" =placebo_date_model_5)

placebo_model_plot <- modelplot(placebo_models_labelled, vcov = vcov_list_placebo_models, size=.1,  coef_omit = c(-2), conf_level=.95, color="black", facet = TRUE, coef_rename = "", ) + geom_vline(xintercept = 0, linetype = 'dashed') + theme_minimal() + xlim(-0.3,0.3) + xlab("Treatment Effect") + theme(axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 8), axis.title.x = element_text(size = 8.5)) + labs(title="Placebo Estimates") + theme(plot.title = element_text(size = 7, hjust = 0.5)) + theme(plot.title = element_text(size = 10))

placebo_model_plot  

combined <- grid.arrange(model_plot, placebo_model_plot, ncol = 2)

# ggsave("Fig6_ITT_Placebo_Estimates.jpeg",combined, width=6 ,height = 4, dpi=600)

#######################################################################################################

#######################################################################################################
#### Study 2: Tables A8 & A9: Full Regression Results, Figure 6 ####

models_coef_map <- c("as.factor(treatment_full)1"="Treatment", "D2pov_age"="Age", "as.factor(D3_education)2"="Secondary, no leaving exam", "as.factor(D3_education)3"= "Secondary, leaving exam", "as.factor(D3_education)4"="University", "c3_religious_attendance"="Rel. Service Attend.", "y2_lib_con"="Lib/Con","as.factor(D1_gender)2"="Female", "r11_income"="Income", "y1_left_right"="Left/Right", "D5_size_settlement"="Size of Settle.", "as.factor(D6_region)2"= "Reg:TN", "as.factor(D6_region)3"= "Reg:TT", "as.factor(D6_region)4"= "Reg:NT", "as.factor(D6_region)5"= "Reg:ZA", "as.factor(D6_region)6"= "Reg:BB", "as.factor(D6_region)7"= "Reg:PR", "as.factor(D6_region)8"= "Reg:KO")

itt_models <- list(outcome_model_bivariate,outcome_model_control_imbalanced,outcome_model_control_all,range2_model,range3_model,range4_model, drop_may3_model, drop_may2_4_model,drop_may1_5_model)

modelsummary(itt_models, output = "latex", vcov = vcov_list_itt_models,  coef_map =models_coef_map , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05,'**'=0.01))


placebo_models <- list(placebo_outcome_model_1,placebo_outcome_model_2,placebo_outcome_model_3,placebo_date_model_1,placebo_date_model_2,placebo_date_model_3, placebo_date_model_4, placebo_date_model_5)


modelsummary(placebo_models, output = "latex", vcov = vcov_list_placebo_models, coef_map =models_coef_map , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05,'**'=0.01))

#######################################################################################################

#######################################################################################################
#### Study 2: Tables A10: Robustness Placebo Models (Including All Controls) ####

data$treatment_full <- as.factor(ifelse(data$datum<"2022-05-03",0,1))

placebo_outcome_all_model_1 <- lm(q10_P_accept_Serb_Romanians ~ as.factor(treatment_full) + D2pov_age + as.factor(D1_gender) + as.factor(D3_education) + r11_income + y1_left_right + y2_lib_con + c3_religious_attendance + as.factor(D6_region) + D5_size_settlement, data=data)
summary(placebo_outcome_all_model_1)

placebo_outcome_all_model_2 <- lm(q10_S_want_anti_corruption~ as.factor(treatment_full) + D2pov_age + as.factor(D1_gender) + as.factor(D3_education) + r11_income + y1_left_right + y2_lib_con + c3_religious_attendance + as.factor(D6_region) + D5_size_settlement, data=data)
summary(placebo_outcome_all_model_2)

placebo_outcome_all_model_3 <- lm(q17_A_beneficial_sending_troops_nato ~ as.factor(treatment_full) + D2pov_age + as.factor(D1_gender) + as.factor(D3_education) + r11_income + y1_left_right + y2_lib_con + c3_religious_attendance + as.factor(D6_region) + D5_size_settlement, data=data)
summary(placebo_outcome_all_model_3)

placebo_date_data <- data[data$datum < as.Date("2022-05-03"), ]

placebo_date_data$treatment_full <- ifelse(placebo_date_data$datum<"2022-04-28",0,1)

placebo_date_all_model_1 <- lm(q16_H_problem_Ukrainian_refugees ~  as.factor(treatment_full) + D2pov_age + as.factor(D1_gender) + as.factor(D3_education) + r11_income + y1_left_right + y2_lib_con + c3_religious_attendance + as.factor(D6_region) + D5_size_settlement, data=placebo_date_data)
summary(placebo_date_all_model_1)

placebo_date_data$treatment_full <- ifelse(placebo_date_data$datum<"2022-04-29",0,1)

placebo_date_all_model_2 <- lm(q16_H_problem_Ukrainian_refugees ~  as.factor(treatment_full) + D2pov_age + as.factor(D1_gender) + as.factor(D3_education) + r11_income + y1_left_right + y2_lib_con + c3_religious_attendance + as.factor(D6_region) + D5_size_settlement, data=placebo_date_data)
summary(placebo_date_all_model_2)

placebo_date_data$treatment_full <- ifelse(placebo_date_data$datum<"2022-04-30",0,1)

placebo_date_all_model_3 <- lm(q16_H_problem_Ukrainian_refugees ~  as.factor(treatment_full) + D2pov_age + as.factor(D1_gender) + as.factor(D3_education) + r11_income + y1_left_right + y2_lib_con + c3_religious_attendance + as.factor(D6_region) + D5_size_settlement, data=placebo_date_data)
summary(placebo_date_all_model_3)

placebo_date_data$treatment_full <- ifelse(placebo_date_data$datum<"2022-05-01",0,1)

placebo_date_all_model_4 <- lm(q16_H_problem_Ukrainian_refugees ~  as.factor(treatment_full) + D2pov_age + as.factor(D1_gender) + as.factor(D3_education) + r11_income + y1_left_right + y2_lib_con + c3_religious_attendance + as.factor(D6_region) + D5_size_settlement, data=placebo_date_data)
summary(placebo_date_all_model_4)

placebo_date_data$treatment_full <- ifelse(placebo_date_data$datum<"2022-05-02",0,1)

placebo_date_all_model_5 <- lm(q16_H_problem_Ukrainian_refugees ~  as.factor(treatment_full) +  D2pov_age + as.factor(D1_gender) + as.factor(D3_education) + r11_income + y1_left_right + y2_lib_con + c3_religious_attendance + as.factor(D6_region) + D5_size_settlement, data=placebo_date_data)
summary(placebo_date_all_model_5)

vcov_list_placebo_all_models <- list(
  vcovHC(placebo_outcome_all_model_1, type = "HC2"),
  vcovHC(placebo_outcome_all_model_2, type = "HC2"),
  vcovHC(placebo_outcome_all_model_3, type = "HC2"),
  vcovHC(placebo_date_all_model_1, type = "HC2"),
  vcovHC(placebo_date_all_model_2, type = "HC2"), 
  vcovHC(placebo_date_all_model_3, type = "HC2"), 
  vcovHC(placebo_date_all_model_4, type = "HC2"), 
  vcovHC(placebo_date_all_model_5, type = "HC2")
)

placebo_models_all <- list(placebo_outcome_all_model_1, placebo_outcome_all_model_2, placebo_outcome_all_model_3, placebo_date_all_model_1, placebo_date_all_model_2, placebo_date_all_model_3, placebo_date_all_model_4, placebo_date_all_model_5)

modelsummary(placebo_models_all, vcov=vcov_list_placebo_all_models, output = "latex", coef_map =models_coef_map , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05,'**'=0.01))

#######################################################################################################

#######################################################################################################
#### Study 2: Figure 7: Conditional Average Treatment Effects: Compliance Moderators ####

data$treatment_full <- as.factor(ifelse(data$datum<"2022-05-03",0,1))

eu_informed_imbalanced_moderation_model <- lm(q16_H_problem_Ukrainian_refugees ~ as.factor(treatment_full)*q18_D_feel_informed_EU + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=data)
summary(eu_informed_imbalanced_moderation_model)

eu_informed_all_moderation_model <- lm(q16_H_problem_Ukrainian_refugees ~ as.factor(treatment_full)*q18_D_feel_informed_EU + D2pov_age + as.factor(D1_gender) + as.factor(D3_education) + c3_religious_attendance + r11_income + y1_left_right + y2_lib_con   + D5_size_settlement, data=data)
summary(eu_informed_all_moderation_model)

slopes(eu_informed_imbalanced_moderation_model, vcov="HC2", variable ="treatment_full",  by="q18_D_feel_informed_EU")

slopes(eu_informed_imbalanced_moderation_model, 
       variable = "treatment_full",
       newdata = datagrid(q18_D_feel_informed_EU = c(7.11)))


base_plot <- plot_slopes(
  model = eu_informed_imbalanced_moderation_model,
  variable = "treatment_full",
  by = "q18_D_feel_informed_EU", 
  vcov="HC2",
  conf_level = 0.95,
  gray = TRUE
)

plot_data <- ggplot_build(base_plot)
x_range_original <- c(1, 10)  
x_range <- c(x_range_original[1] - 0.25, x_range_original[2] + 0.25)
x_breaks <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)  

main_plot <- plot_slopes(
  model = eu_informed_imbalanced_moderation_model,
  variable = "treatment_full",
  by = "q18_D_feel_informed_EU",  
  vcov = "HC2", 
  conf_level = 0.95,
  gray = TRUE
) +
  scale_x_continuous(breaks = x_breaks, limits = x_range) + 
  scale_y_continuous(limits = c(-0.65, 0.05), breaks = seq(-0.6, 0.0, by = 0.1)) + 
  theme_minimal() +
  xlab("") +  # Remove x-axis label from main plot
  ylab("Marginal Effect of Treatment") +
  labs(title = "") +
  theme(
    axis.text.x = element_blank(),  
    axis.text.y = element_text(size = 6.5),
    axis.title.x = element_blank(),  
    axis.title.y = element_text(size = 8),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.title = element_blank(),
    plot.margin = margin(5, 5, 0, 5)  # Remove bottom margin
  )

hist_plot_grouped <- ggplot(data, aes(x = q18_D_feel_informed_EU, fill = as.factor(treatment_full))) +
  geom_histogram(alpha = 0.7, bins = 20, position = position_dodge(width = .45)) +  # Stagger the bars
  scale_x_continuous(breaks = x_breaks, limits = x_range) + 
  scale_y_continuous(limits = c(0, 625), breaks = seq(0, 625, by = 200)) + 
  scale_fill_manual(
    values = c("0" = "black", "1" = "gray40"),
    labels = c("0" = "Control", "1" = "Treatment"),
    name = ""
  ) +
  theme_minimal() +
  xlab("Feel Informed About EU") +
  ylab("Count") +
  theme(
    axis.text.x = element_text(size = 6),
    axis.text.y = element_text(size = 5),
    axis.title.x = element_text(size = 7),
    axis.title.y = element_text(size = 7),
    panel.grid.minor.y = element_blank(), 
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.text = element_text(size = 6),
    legend.key.size = unit(0.3, "cm"),
    plot.margin = margin(0, 5, 5, 5)
  )

marginal_effect_plot_eu_informed <- main_plot / hist_plot_grouped +
  plot_layout(heights = c(3, 1))

marginal_effect_plot_eu_informed


media_trust_imbalanced_moderation_model <- lm(q16_H_problem_Ukrainian_refugees ~ as.factor(treatment_full)*q4_R_trust_standard_media + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=data)
summary(media_trust_imbalanced_moderation_model)


media_trust_all_moderation_model <- lm(q16_H_problem_Ukrainian_refugees ~ as.factor(treatment_full)*q4_R_trust_standard_media + D2pov_age + as.factor(D1_gender) + as.factor(D3_education) + c3_religious_attendance  + r11_income + y1_left_right + y2_lib_con   + D5_size_settlement, data=data)
summary(media_trust_all_moderation_model)

avg_slopes(media_trust_imbalanced_moderation_model, vcov="HC2", variable ="treatment_full",  by="q4_R_trust_standard_media")

slopes(media_trust_imbalanced_moderation_model, 
       variable = "treatment_full",
       newdata = datagrid(q4_R_trust_standard_media = c(1, 4)))


base_plot <- plot_slopes(
  model = media_trust_imbalanced_moderation_model,
  variable = "treatment_full",
  by = "q4_R_trust_standard_media",
  vcov = "HC2", 
  conf_level = 0.95,
  gray = TRUE
)

plot_data <- ggplot_build(base_plot)
x_range_original <- c(1, 4)  
x_range <- c(x_range_original[1] - 0.15, x_range_original[2] + 0.15)
x_breaks <- plot_data$layout$panel_scales_x[[1]]$breaks

main_plot <- plot_slopes(
  model = media_trust_imbalanced_moderation_model,
  variable = "treatment_full",
  by = "q4_R_trust_standard_media",   
  vcov = "HC2", 
  conf_level = 0.95,
  gray = TRUE
) +
  scale_x_continuous(breaks = x_breaks, limits = x_range) + 
  scale_y_continuous(limits = c(-0.65, 0.05), breaks = seq(-0.6, 0.0, by = 0.1)) + 
  theme_minimal() +
  xlab("") +  
  ylab("") +
  labs(title = "") +
  theme(
    axis.text.x = element_blank(),  
    axis.text.y = element_text(size = 6.5),
    axis.title.x = element_blank(),  
    axis.title.y = element_text(size = 8),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.title = element_blank(),
    plot.margin = margin(5, 5, 0, 5)  
  )


hist_plot_grouped <- ggplot(data, aes(x = q4_R_trust_standard_media, fill = as.factor(treatment_full))) +
  geom_histogram(alpha = 0.7, bins = 20, position = position_dodge(width = 0.175)) +  # Stagger the bars
  scale_x_continuous(breaks = x_breaks, limits = x_range) + 
  scale_y_continuous(limits = c(0, 625), breaks = seq(0, 625, by = 200)) + 
  scale_fill_manual(
    values = c("0" = "black", "1" = "gray40"),
    labels = c("0" = "Control", "1" = "Treatment"),
    name = ""
  ) +
  theme_minimal() +
  xlab("Media Trust") +
  ylab("") +
  theme(
    axis.text.x = element_text(size = 6),
    axis.text.y = element_text(size = 5),
    axis.title.x = element_text(size = 7),
    axis.title.y = element_text(size = 7),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.text = element_text(size = 6),
    legend.key.size = unit(0.3, "cm"),
    plot.margin = margin(0, 5, 5, 5)
  )

marginal_effect_media_trust <- main_plot / hist_plot_grouped +
  plot_layout(heights = c(3, 1))

marginal_effect_media_trust

combined_marginal_effects <- ggarrange(marginal_effect_plot_eu_informed,marginal_effect_media_trust, ncol=2, common.legend = TRUE, legend="bottom")

combined_marginal_effects

# ggsave("Fig7_Compliance_Moderators.jpeg", plot = combined_marginal_effects, width=6 ,height = 5, dpi=600)

#######################################################################################################

#######################################################################################################
#### Study 2: Table A11: Full Regression Results and Robustness Compliance Models Figure 7 Main Text ####


itt_moderation_models <- list(eu_informed_imbalanced_moderation_model,media_trust_imbalanced_moderation_model, eu_informed_all_moderation_model, media_trust_all_moderation_model)

itt_moderation_coef_map <- c("as.factor(treatment_full)1"="Treatment", "q18_D_feel_informed_EU"="Informed EU", "q4_R_trust_standard_media"="Media Trust", "as.factor(treatment_full)1:q18_D_feel_informed_EU" = "Treatment*Informed EU", "as.factor(treatment_full)1:q4_R_trust_standard_media"="Treatment*Media Trust", "as.factor(treatment_full)1:y2_lib_con"="Treatment*Lib/Con", "D2pov_age"="Age", "as.factor(D3_education)2"="Secondary, no leaving exam", "as.factor(D3_education)3"= "Secondary, leaving exam", "as.factor(D3_education)4"="University", "c3_religious_attendance"="Rel. Service Attend.", "y2_lib_con"="Lib/Con","as.factor(D1_gender)2"="Female", "r11_income"="Income", "y1_left_right"="Left/Right", "D5_size_settlement"="Size of Settle.", "as.factor(D6_region)2"= "Reg:TN", "as.factor(D6_region)3"= "Reg:TT", "as.factor(D6_region)4"= "Reg:NT", "as.factor(D6_region)5"= "Reg:ZA", "as.factor(D6_region)6"= "Reg:BB", "as.factor(D6_region)7"= "Reg:PR", "as.factor(D6_region)8"= "Reg:KO")

modelsummary(itt_moderation_models, output = "latex", vcov = "HC2", coef_map =itt_moderation_coef_map , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*'=0.05,'**' = .01))

#######################################################################################################

#######################################################################################################
#### Study 2: Figure 8: Conditional Average Treatment Effects: Ideology and Locus of Control ####

lib_con_imbalanced_moderation_model <- lm(q16_H_problem_Ukrainian_refugees ~ as.factor(treatment_full)*y2_lib_con + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=data)
summary(lib_con_imbalanced_moderation_model)

avg_slopes(lib_con_imbalanced_moderation_model, variable= "treatment_full", by="y2_lib_con", vcov="HC2")

sd(data$q16_H_problem_Ukrainian_refugees)

lib_con_all_moderation_model <- lm(q16_H_problem_Ukrainian_refugees ~ as.factor(treatment_full)*y2_lib_con + D2pov_age + as.factor(D1_gender) + as.factor(D3_education) + c3_religious_attendance + r11_income + y1_left_right + y2_lib_con + D5_size_settlement, data=data)
summary(lib_con_all_moderation_model)

base_plot <- plot_slopes(
  model = lib_con_imbalanced_moderation_model,
  variable = "treatment_full",
  by = "y2_lib_con",    
  vcov = "HC2", 
  conf_level = 0.95,
  gray = TRUE
)

plot_data <- ggplot_build(base_plot)
x_range_original <- c(1, 5) 
x_range <- c(x_range_original[1] - 0.25, x_range_original[2] + 0.25)
x_breaks <- plot_data$layout$panel_scales_x[[1]]$breaks

main_plot <- plot_slopes(
  model = lib_con_imbalanced_moderation_model,
  variable = "treatment_full",
  by = "y2_lib_con",    
  vcov = "HC2", 
  conf_level = 0.95,
  gray = TRUE
) +
  scale_x_continuous(breaks = x_breaks, limits = x_range) +  
  scale_y_continuous(limits=c(-0.6,0.2), breaks=seq(-.6, 0.2, by = 0.2), labels = label_number(accuracy = 0.1)) +
  theme_minimal() +
  xlab("") +  # Remove x-axis label from main plot
  ylab("Marginal Effect of Treatment") +
  labs(title = "") +
  theme(
    axis.text.x = element_blank(),  
    axis.text.y = element_text(size = 6.5),
    axis.title.x = element_blank(),  
    axis.title.y = element_text(size = 8),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.title = element_blank(),
    plot.margin = margin(5, 5, 0, 5)  # Remove bottom margin
  )

hist_plot_grouped <- ggplot(data, aes(x = y2_lib_con, fill = as.factor(treatment_full))) +
  geom_histogram(alpha = 0.7, bins = 20, position = position_dodge(width = 0.22)) +  # Stagger the bars
  scale_x_continuous(breaks = x_breaks, limits = x_range) + 
  scale_y_continuous(limits = c(0, 900), breaks = seq(0, 900, by = 200)) +
  scale_fill_manual(
    values = c("0" = "black", "1" = "gray40"),
    labels = c("0" = "Control", "1" = "Treatment"),
    name = ""
  ) +
  theme_minimal() +
  xlab("Lib/Con") +
  ylab("Count") +
  theme(
    axis.text.x = element_text(size = 6),
    axis.text.y = element_text(size = 5),
    axis.title.x = element_text(size = 7),
    axis.title.y = element_text(size = 7),
    panel.grid.minor.y = element_blank(), 
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.text = element_text(size = 6),
    legend.key.size = unit(0.3, "cm"),
    plot.margin = margin(0, 5, 5, 5)
  )

marginal_effect_plot_lib_con_grouped <- main_plot / hist_plot_grouped +
  plot_layout(heights = c(3, 1))

marginal_effect_plot_lib_con_grouped


locus_imbalanced_moderation_model <- lm(q16_H_problem_Ukrainian_refugees ~ as.factor(treatment_full)*q10_B_who_power_does_not_matter + D2pov_age + as.factor(D3_education) + c3_religious_attendance + y2_lib_con, data=data)
summary(locus_imbalanced_moderation_model)

avg_slopes(locus_imbalanced_moderation_model, variable= "treatment_full", by="q10_B_who_power_does_not_matter", vcov="HC2")


locus_all_moderation_model <- lm(q16_H_problem_Ukrainian_refugees ~ as.factor(treatment_full)*q10_B_who_power_does_not_matter + D2pov_age + as.factor(D1_gender) + as.factor(D3_education) + c3_religious_attendance + r11_income + y1_left_right + y2_lib_con + D5_size_settlement, data=data)
summary(locus_all_moderation_model)


base_plot <- plot_slopes(
  model = locus_imbalanced_moderation_model,
  variable = "treatment_full",
  by = "q10_B_who_power_does_not_matter",    
  vcov = "HC2", 
  conf_level = 0.95,
  gray = TRUE
)

plot_data <- ggplot_build(base_plot)
x_range_original <- c(1, 4) 
x_range <- c(x_range_original[1] - 0.25, x_range_original[2] + 0.25)
x_breaks <- plot_data$layout$panel_scales_x[[1]]$breaks

main_plot <- plot_slopes(
  model = locus_imbalanced_moderation_model,
  variable = "treatment_full",
  by = "q10_B_who_power_does_not_matter",    
  vcov = "HC2", 
  conf_level = 0.95,
  gray = TRUE
) +
  scale_x_continuous(breaks = x_breaks, limits = x_range) +  
  scale_y_continuous(limits=c(-0.6,0.2), breaks=seq(-.6, 0.2, by = 0.2), labels = label_number(accuracy = 0.1)) +
  theme_minimal() +
  xlab("") +  # Remove x-axis label from main plot
  ylab("") +
  labs(title = "") +
  theme(
    axis.text.x = element_blank(),  
    axis.text.y = element_text(size = 6.5),
    axis.title.x = element_blank(),  
    axis.title.y = element_text(size = 8),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.title = element_blank(),
    plot.margin = margin(5, 5, 0, 5)  # Remove bottom margin
  )

hist_plot_grouped <- ggplot(data, aes(x = q10_B_who_power_does_not_matter, fill = as.factor(treatment_full))) +
  geom_histogram(alpha = 0.7, bins = 20, position = position_dodge(width = 0.175)) +  # Stagger the bars
  scale_x_continuous(breaks = x_breaks, limits = x_range) + 
  scale_y_continuous(limits = c(0, 900), breaks = seq(0, 900, by = 200)) +
  scale_fill_manual(
    values = c("0" = "black", "1" = "gray40"),
    labels = c("0" = "Control", "1" = "Treatment"),
    name = ""
  ) +
  theme_minimal() +
  xlab("State Locus of Control") +
  ylab("") +
  theme(
    axis.text.x = element_text(size = 6),
    axis.text.y = element_text(size = 5),
    axis.title.x = element_text(size = 7),
    axis.title.y = element_text(size = 7),
    panel.grid.minor.y = element_blank(), 
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.text = element_text(size = 6),
    legend.key.size = unit(0.3, "cm"),
    plot.margin = margin(0, 5, 5, 5)
  )

marginal_effect_plot_locus_grouped <- main_plot / hist_plot_grouped +
  plot_layout(heights = c(3, 1))

marginal_effect_plot_locus_grouped


combined_marginal_effects_2 <- ggarrange(marginal_effect_plot_lib_con_grouped,marginal_effect_plot_locus_grouped, ncol=2, common.legend = TRUE, legend="bottom")

combined_marginal_effects_2

# ggsave("Fig8_LoC_Moderators.jpeg", plot = combined_marginal_effects_2, width=6 ,height = 5, dpi=600)

#######################################################################################################

#######################################################################################################
#### Study 2: Table A12: Full Regression Results CATE Models Figure 8 Main Text ####

moderation_models <- list(lib_con_imbalanced_moderation_model,lib_con_all_moderation_model, locus_imbalanced_moderation_model, locus_all_moderation_model)

cm_moderation_models <- c("as.factor(treatment_full)1"="Treatment", "y2_lib_con"="Lib/Con", "q10_B_who_power_does_not_matter"="State LoC", "as.factor(treatment_full)1:y2_lib_con"="Treatment*Lib/Con", "as.factor(treatment_full)1:q10_B_who_power_does_not_matter" = "Treatment*State LoC", "D2pov_age"="Age", "as.factor(D3_education)2"="Secondary, no leaving exam", "as.factor(D3_education)3"= "Secondary, leaving exam", "as.factor(D3_education)4"="University", "c3_religious_attendance"="Rel. Service Attend.","as.factor(D1_gender)2"="Female", "r11_income"="Income", "y1_left_right"="Left/Right", "D5_size_settlement"="Size of Settle.", "as.factor(D6_region)2"= "Reg:TN", "as.factor(D6_region)3"= "Reg:TT", "as.factor(D6_region)4"= "Reg:NT", "as.factor(D6_region)5"= "Reg:ZA", "as.factor(D6_region)6"= "Reg:BB", "as.factor(D6_region)7"= "Reg:PR", "as.factor(D6_region)8"= "Reg:KO")

modelsummary(moderation_models, vcov = "HC2", output = "latex", coef_map =cm_moderation_models , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*'=0.05,'**' = .01))

#######################################################################################################
